Matter power spectrum from a Lagrangian-space regularization of perturbation theory 
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We present a new approach to compute the matter density power spectrum, from large linear 
scales to small highly nonlinear scales. Instead of explicitly computing a partial series of high- 
order diagrams, as in perturbative resummation schemes, we embed the standard perturbation 
theory within a realistic nonlinear Lagrangian-space ansatz. We also point out that an "adhesion- 
like" regularization of the shell crossing regime is more realistic than a "Zel'dovich-like" behavior, 
where particles freely escape to infinity. This provides a "cosmic web" power spectrum with good 
small-scale properties that provides a good matching with a halo model on mildly nonlinear scales. 
We obtain a good agreement with numerical simulations on large scales, better than 3% for k < 
l^Mpc"^, and on small scales, better than 10% for k < lO/iMpc"^, a.t z > 0.35, which improves 
over previous methods. 

PACS numbers: 98.80.-k, 98.65.Dx 



I. INTRODUCTION 

The description of gravitational clustering can follow 
two different frameworks, the Eulerian approach, where 
one studies density and velocity fields on a fixed grid, 
and the Lagrangian approach, where one studies particle 
trajectories themselves [l|. 

The Eulerian approach is more convenient for some 
practical purposes because we observe density fields and 
not trajectories (and it is also better suited to bary- 
onic physics, which involves pressure and temperature 
fields). However, it has the theoretical disadvantage that 
one should in principle work with the phase-space distri- 
bution function /(x, v;t), which involves seven variables 
and is very heavy for numerical and analytical compu- 
tations. Then, in practice most analytical approaches 
are based on the fluid approximation, which replaces the 
Vlasov equation by hydrodynamical equations for the 
density and velocity fields p(x, t) and v(x, t), which are 
much easier to handle. However, this is only exact in 
the single-stream regime and these equations of motion 
themselves break down after shell crossing. This makes 
the matching between the perturbative and nonperturba- 
tive regimes rather difficult. More precisely, because the 
dynamics is not well defined (these equations of motion 
are not sufficient to determine the evolution after shell 
crossing) the asymptotic regime at high wavenumber of 



the power spectrum obtained within this framework is 
not obvious and not clearly related to physical behav- 
iors. 

In contrast, in the Lagrangian approach the trajecto- 
ries x(q, t), where particles are identified by their initial 
position q, are always well defined and the behavior after 
shell crossing is (at least implicitly) defined within each 
peculiar scheme. This offers the prospect of a more con- 
venient matching to the highly nonlinear regime. More- 
over, in contrast to the Eulerian case the Lagrangian ap- 
proach is not sensitive to the "sweeping effect" associated 
with almost uniform translations by long wavelengths of 
the velocity field 0411 (i-e., this can be automatically 
removed as a random uniform shift) and it is a more di- 
rect probe of the structures of the density field. From an 
observational point of view, this framework is also more 
convenient to include redshift-space distortions. How- 
ever, to compare theoretical predictions with observa- 
tions one typically needs to compute again density fields 
and the derivation of the latter from the particle trajec- 
tories usually introduces additional approximations and 
problems. 

This latter disadvantage explains why most perturba- 
tive schemes that have been recently developed follow the 
Eulerian framework [Ml. They provide a good match 
with numerical simulations on quasilinear scales by using 
resummation schemes that are exact up to one or two- 
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loop order and improve over the standard perturbative 
approach by including partial resummations of higher- 
order diagrams. However, it has been difficult so far to 
obtain a good model up to the highly nonlinear scales 
(by embedding the perturbative predictions within a phe- 
nomenological halo model) as one usually underestimates 
the power spectrum on transition scales [ll|, [l^. This 
discrepancy is even worse for Lagrangian schemes, which 
usually give rise to a perturbative power spectrum that 
decays faster than the linear power spectrum at high k 
[Til IigI . [TtI . whereas the Eulerian schemes can show a 
variety of behaviors, a similar fast decay, a convergence 
back to the linear power, or a higher tail. 

In this paper, we follow two main motivations: 

- develop a Lagrangian-based approach, 

- obtain a good matching with the highly nonlinear 
regime. 

Indeed, a Lagrangian-based approach is expected to 
provide a better starting point for eventually tackling 
redshift-space distortions and biasing (although we do 
not address these points here). Moreover, it is not sen- 
sitive to the "sweeping effect" , so that the modeling can 
be directly focused on density structures, that is, rela- 
tive displacements. Next, because our final goal is to 
obtain a unified description from linear to highly nonlin- 
ear scales, it is necessary to pay attention to transition 
scales. Moreover, the requirement of a good matching to 
nonlinear scales is likely to shed light on the asymptotic 
behavior that should be satisfied by good perturbative 
schemes. 

To reach these goals, we advocate a new perspective 
while drawing on some previous works. Our main new 
ideas are the following: 

a) instead of building a perturbative resummation 
scheme, which explicitly computes a partial series of high- 
order diagrams, we build a self-consistent ansatz that we 
next make consistent with standard perturbation theory 
up to some order (one- loop order in this paper), 

b) instead of the usual "Zel'dovich-like" behavior in 
the shell-crossing regime we impose an "adhesion-like" 
behavior. 

These two strategies are motivated by the wish to ob- 
tain a good behavior up to nonlinear scales. 

First, most perturbative schemes are not guaranteed to 
provide at each order self-consistent predictions, whether 
they involve a sharp truncation at a finite order (as in 
standard perturbation theory) or include partial resum- 
mations of high-order terms. In other words, they do 
not guarantee that their predictions can be realized by a 
physical density field (such that the density is real and 
positive) and one may encounter unphysical behaviors 
at high k (such as a negative power spectrum). This is 
not necessarily a problem, if one restricts to the range 
of validity of the perturbative scheme. However, since 
we intend to build a unified model that applies up to 
nonlinear scales, we need to obtain as a building block a 
perturbative power spectrum that remains well behaved 
in all regimes. A second reason is that physical self- 



consistency (at least to some degree) can be expected 
to ensure good convergence properties or at least a rea- 
sonable high-fc tail. For instance, this clearly rules out 
the bad behavior of the standard Eulerian perturbation 
theory, where higher-order contributions grow increas- 
ingly fast on small scales and lead to power spectra that 
can change sign with the truncation order. Then, we 
can hope that a better-controlled high-A: behavior and a 
better convergence can in turn improve the accuracy on 
quasilinear scales. 

The method that we develop to achieve this goal is 
to first build an ansatz for the power spectrum that is 
always well behaved and next to tune its parameters (in 
our case, the skewness of the longitudinal displacement) 
so that its perturbative expansion matches the standard 
perturbative expansion up to the required order (here up 
to second order over Pl)- 

Second, the "Zel'dovich-like" behavior of usual La- 
grangian perturbative schemes, where particles escape to 
infinity after shell crossing, leads to a fast decay of the 
power spectrum on nonlinear scales. This is due to the 
erasing of intermediate-scale structures such as pancakes 
soon after their formation, because particles do not re- 
main trapped within potential wells. As recalled above, 
this may be the source of a significant underestimation of 
the power spectrum on transition scales. The "adhesion 
model" introduced in [l^, where particles stick together 
after shell crossing, seems a more realistic approxima- 
tion. In particular, this provides a good description of 
the cosmic web. Because it is not easy to implement 
the exact adhesion model within a model for the power 
spectrum we use a simplified ansatz, inspired from [l9j . 
to include some form of "sticking" of particle pairs after 
shell crossing. This models the formation of pancakes 
and yields additional power at high k, which should be 
more realistic and more accurate than previous methods. 

These ingredients allow us to build a model for the 
"cosmic web" power spectrum, associated with large and 
intermediate scale structures (bulk fiows, voids, and pan- 
cackes). Next, following [UlliBl, we combine these results 
with a phenomenological halo model to extend our model 
up to highly nonlinear scales, associated with inner halo 
regions. 

This paper is organized as follows. After recalling the 
expression of the matter density power spectrum in a 
Lagrangian-space framework in Sec. [Hi we describe in 
Sec. mil our model for the "cosmic web" power spec- 
trum. This corresponds to both the large-scale perturba- 
tive regime and the nonperturbative intermediate-scale 
regime associated with pancakes. Next, we explain in 
Sec. lIVI how we combine this cosmic web power spectrum 
with a halo model to extend our model to highly nonlin- 
ear scales, associated with inner halo regions. Then, we 
compare our model with numerical simulations in Sec. IVl 
and wc conclude in Sec. IVII 

The reader who is only interested in our final numerical 
results may directly go to Sec. |Vl 

In this paper, we mainly focus on a fiat ACDM cosmol- 
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ogy derived from the five year observation by the WMAP 
satelhte [l^l ■ We measure some of the ingredients in our 
model from five out of 60 reaUzations of N-body simula- 
tions performed in adopting this cosmological model. 
The numerical data for the power spectrum at large scale 
(fc < Ih Mpc^^) is taken from the full 60 simulations in 
that paper, while we refer to a higher-resolution simu- 
lation by for smaller scales. Wc finally show how 
well our model can reproduce the power spectrum from 
N-body simulations done in two other cosmologies. 



II. MATTER DENSITY POWER SPECTRUM IN 
A LAGRANGIAN FRAMEWORK 

In a Lagrangian framework one considers the trajec- 
tories x(q, t) of all particles, of initial Lagrangian co- 
ordinates q and Eulerian coordinates x at time t. In 
particular, at any given time t, this defines a mapping, 
q I— > X, from Lagrangian to Eulerian space, which fully 
determines the Eulerian density field /9(x) through the 
conservation of matter. 



p(x) dx = pdq. 



(1) 



where p is the mean comoving matter density of the Uni- 
verse and we work in comoving coordinates. Then, defin- 
ing the density contrast as 



<5(x,t) 



p(x,t) - p 



and its Fourier transform as 

dx 



<5(k) 



<5(x). 



(2) 



(3) 



,-ikq 



(4) 



(27r)3 

one obtains from Eq.(IT]) 

*«-/(^('-"* 

Next, defining the density power spectrum as 

(^(ki)5(k2)) = Soik, + k2)P(fcl), (5) 

we obtain from Eq. ^ , using statistical homogeneity [2l], 



(6) 



where we introduced the Lagrangian-space and Eulerian- 
space separations Aq and Ax, 



Aq = q2 - qi , Ax = X2 - xi , 



(7) 



between two particles qi and q2. The expression ([6]) is 
fully general since it is a simple consequence of the matter 
conservation Eq.(IT]), and of statistical homogeneity. In 
particular, it holds for any dynamics, such as the one 
associated with the Zel'dovich approximation psj . where 
the mapping x(q) is given by the linear displacement 
field. 



III. LARGE-SCALE POWER SPECTRUM 
(COSMIC WEB) 

Our final goal is to build an unified model for the power 
spectrum that applies from linear to highly nonlinear 
scales, which requires at some level a phenomenological 
model to describe small scales. As described in Sec. IIVI 
below, following [ll|, [l5| we use a halo model to combine 
large-scale perturbative schemes with phenomenological 
approaches. However, wc first focus on our perturba- 
tive approach, where we neglect the formation of virial- 
ized halos. Thus, we consider in this section the power 
spectrum associated with large-scale bulk flows, that is, 
the formation of the cosmic web, disregarding small scale 
structures such as inner halo regions. This actually in- 
cludes two regime: i) the very large scales where shell 
crossing can be neglected, and ii) the intermediate scales, 
associated with pancakes or filaments, where shell cross- 
ing comes into play to shape the cosmic web but where 
the internal structure of pancakes and filaments can be 
neglected. The first regime can be described by pertur- 
bation theory and it is considered in Sees . IIII XI and HIl B I 
Next, we tackle the second regime in Sec. IIIICI 



A. Zel'dovich power spectrum 

If all particles pairs can be described by the same 
perturbative framework, the power spectrum ([5]) can be 
written as 



P{k) = 



dAq 

dAq 

(27r)3 



(e 



ikAx\ 



exp 



^ ((ik ■ Ax)"), 

n=l 



n 



(8) 
(9) 



In the first line we used the fact that the integral of the 
last term in Eq.® gives a Dirac factor Snik.) that can 
be discarded for fc > and in the second line we used the 
usual expansion over cumulants. 

In the well-known Zel'dovich approximation (23| we 
use the linear prediction for the Eulerian separation Ax. 
Then, for Gaussian initial conditions only the first and 
second-order cumulants are nonzero. Introducing the dis- 
placement field x(q, t) = q -I- *(q, t), at linear order 
wc have the longitudinal and transverse variances (with 
respect to the direction defined by the initial Lagrangian 
separation Aq) 

af^{Aq) = 2 J dk[l-cos(fciA(7)]^Pi.(fc), (10) 
^i(Ag) = 2 J dk[l-cos(A:iAg)]^p^(fc). (u) 

Here we took Aq = {Aq)ei along the first axis, ei, and 
wc labeled the two transverse axis as e2 and e^. Thus, 
aj = {(A^m?) = {(A^Li?) and a'i = {iA^L2)'') = 
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FIG. 1: Matter density power spectrum a.t z = 0.35 (mul- 
tiplied by a factor (kh)^'^ to decrease the vertical range). 
We show the results from N-body simulations (data points), 
the linear power spectrum Pl, the standard one-loop pre- 



diction Piioop, the one- loop prediction Pnoop lfT8|) within the 
Zel'dovich approximation, and the full nonlinear Zel'dovich 
power spectrum (|16|) . 



{{A'i> Ls)"^) . Introducing the integrals 



Iiiq)^Y Jo '^'^ P^ik) Mqk), 



(12) 



and the variance of the linear one-point displacement 
along one dimension, 



a2 = i(|*^(q)p)=/o(0), 



(13) 



we have: 



af(Ag) = 2al~2h{Aq)+Al2(Aq), (14) 
aliAq) = 2al-2Io{Aq)-2l2{Aq). (15) 

Then, in the Zel'dovich approximation, where x = q + 
*L, Eq.(ini) reads as 



(27r)3 



, (16) 



where /i = (k • Aq)/ (kAq). (By symmetry, we always 
have (A\I/) = 0.) At large distance, the linear displace- 
ments '4'(qi) and 4'(q2) of the two particles become in- 
dependent and crjj and a'j_ converge to 2ct^ , in agreement 
with Eqs.dHIl-dini). Then, Eq.idl]) is formally divergent 
at large q, as it contains a Dirac factor ^^(k), and for 
numerical computations it is convenient to compute in a 
separate fashion low-order terms. 

By expanding the exponential ^ over powers of Pl 
one recovers the standard Eulerian perturbation theory. 



Within the Zel'dovich approximation, we obtain for in- 
stance up to second order [23| 



p2(fc)=Pi(fc) + pZ„,p(fc) 



(17) 



with 



Pi 



iioop(fc) = -k'alPLik) + j dkidka fe(ki + ka - k 
(k.ki)2(k.k2)2 



-PLikl)PL{k2). 



(18) 



This is different from the exact one-loop contribution 
Piioop generated by the actual gravitational dynamics 
because at this order we should include the third-order 
cumulants and loop corrections to the second-order cu- 
mulants (the cumulant of order n scales as ((A'S')")^ oc 
P2^^ at leading order). 

We compare in Fig. [T] the power spectra obtained 
from the true gravitational dynamics (at linear and one- 
loop order, as well as the full nonlinear power spectrum 
from N-body simulations) with those associated with the 
Zel'dovich approximation (at one-loop order and for the 
full nonlinear power spectrum). As is well known, the 
standard Eulerian one-loop prediction Piioop improves 
the agreement with simulations as compared with lin- 
ear theory on large scales, but quickly gives too much 
power and is badly behaved at high k. In contrast, the 
nonlinear Zel'dovich power spectrum P^ decays faster 
than the linear power at high k [2l], |2^. For instance, 
for a power-law linear power spectrum Phik) oc fc" with 
-3 < n < -1 we have P'^{k) - fc-3+3(n+3)/(n-n) j^jgj^ 

k [3]. This is due to the fact that within the approxi- 
mation of linear trajectories particles keep moving along 
straight lines after shell crossing so that pancakes and fil- 
aments quickly fatten and dissolve after their formation. 
This also leads to a one-loop Zel'dovich prediction 
that is smaller than the linear power spectrum on weakly 
nonlinear scales, where it is meaningful and follows the 
full nonlinear power P^. Therefore, at one- loop order 
the Zel'dovich approximation is actually worse than the 
simple linear approximation. 

The figure[T]teaches us two things if we wish to describe 
the power spectrum in a Lagrangian framework: 

i) we must explicitly modify the Zel'dovich approxima- 
tion at one-loop order to (partly) cure the underestima- 
tion of the power spectrum on weakly nonlinear scales, 

ii) we must modify the behavior after shell crossing to 
cure the high-Zc damping. 



B. Perturbative Lagrangian ansatz 

1. Factorized ansatz 



Before tackling the high-fc regime in section IIII CI 
(point ii) above) , we address the weakly nonlinear regime 
in this section (point i)), where we disregard shell cross- 
ing. Thus, we look for a Lagrangian ansatz that improves 



5 



over (|16p on weakly nonlinear scales and shows a pertur- 
bative expansion over integer powers of Pi^ as in the true 
gravitational dynamics. 

By symmetry the means ((A^'n )^A'>I'_l)c and 
((A*_l)'^)c are zero and at order P| we are left 
with the two third-order cumulants ((A5'j|)'^)c and 
(A^*!! (A5'j^)^)c. In particular, the latter average shows 
that beyond the Gaussian order (fTB|) the longitudinal 
and transverse displacements are generically correlated. 
Nevertheless, because our goal is to build a simple 
ansatz for Eq.Q we consider a simplified model where 
the longitudinal and transverse displacements are uncor- 
related. Then, since (A4'^)c is of order P| we keep the 
linear Gaussian for the transverse part and we look for 
an ansatz of the form 



P\k) 



dAq 

(2^ 



(19) 



where Aa;|| = Aq + is the longitudinal Eulerian sep- 
aration and we need to specify a model for the mean 
Thus, the difference from the Zel'dovich approxi- 
mation (|16p is that we go beyond the Gaussian for the 
longitudinal part. Apart from simplicity and ordering in 
perturbation theory, the reason why we focus on the lon- 
gitudinal part is that it should be more sensitive than 
the transverse part to the progress of nonlinear cluster- 
ing, even at a qualitative level. Indeed, by symmetry 
the probability distribution of the relative displacement 
along any transverse direction, 7'(A\E'j^), remains even 
and peaks at A^I*^ = 0, so that on a qualitative level 
its shape remains similar to a Gaussian. In contrast, the 
probability distribution of the longitudinal Eulerian sep- 
aration, P(Ax||), which in the linear regime is a Gaussian 
centered on Aa;|| — Aq, develops an asymmetry (skew- 
ness) that is characteristic of the gravitational dynamics. 
Moreover, following [l^, negative values of Ax|| are ex- 
pected to be closely associated with shell crossing, so 
that keeping track of Axu will be useful to handle shell- 
crossing effects as explained in Sec. IIII Cl below. 



2. One-loop order 

It is convenient to define the dimensionless Eulerian 
longitudinal separation kii by 



Aa 



= 1 



A*ii 



Ag ^ ' Ag ' 
and its linear variance and its skewncss by 



(20) 




Aq [/7"'Mpc] 



FIG. 2: The linear variances a\\ and ax (in units of /i^^Mpc), 
and (Tk, at z = 0.35. The lines are the linear-theory predic- 
tions (|10|l - (|lip and the points are the results from N-body 
simulations. 



linear theory gives a good estimate of these variances. On 
smaller scales, nonlinearities and nonperturbative effects 
(virialized motions in the shell-crossing regime), which 
we do not consider in this section, increase the ampli- 
tude of these rms displacements. As noticed above, at 
large distance a\\ and a±^ converge to %/2ct„ because the 
motions of the two particles become independent. This 
implies that cr„ decreases as l/(Ag) on large scales. On 
small separations, for CDM power spectra that decrease 
faster than at high k, a\\ and (t_l decrease as Ag and 
(Tk converges to a constant. However, we can see that 
these asymptotic regimes are only reached on very large 
or very small scales. The value of CTk is a better mea- 
sure of nonlinearity than the rms displacements a\\ and 
(Tj^, and (Tk ^ 1 marks the transition to the nonlinear 
regime. From the definition (|20p . in a spherical dynam- 
ics K|| (1 + 6)^^/^ , and deviations of order unity of K|| 
from its mean are associated with density contrasts of 
order unity. In particular, K|| — 1 = —1 corresponds to 
shell crossing and infinite density. In CDM cosmologies 
smaller scales turn nonlinear first and we can check in 
Fig. [2] that (Tk increases on smaller scales. 

To further simplify our ansatz, we keep the linear the- 
ory prediction for the variance (k|)c. Then, the only 
difference at order P| of the power spectrum ([TO]) from 
the Zel'dovich approximation p6p arises from the skew- 
ness Sr, " : 



'II 



(Ag)2 



,2\2 



(A<z) 



((A^||)^)c 
((A*||)')? 



(21) 

We show the linear variances a\\ and (Tk|| (as well as 
the transverse counterparts cr^ and ir^^ = (T^/{/S.q)) in 
Fig. m We can check that above 3/i~^Mpc, at z = 0.35, 



with 



PLo^{k) = Pl,,^{k) + PsAk), 



(22) 



PsM = / ^ sin(fc/.Ag) (fc/.Ag)3 ^i^. (23) 
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The integration over angles gives 



dAq 



(Aq)^ 



27r2 

X cos{kAq) + [{kAq)- 



[6 - (kAq)' 



6 



, sin(A:Ao) 1 , , 
2] 34^}. (24) 



Then, to complete our model at one-loop order we must 
set the skewness " . A first choice would be to use the 
prediction from lowest-order perturbation theory. Then, 
going up to second order in Lagrangian perturbation the- 
ory ("2LPT", see H-Hl), we obtain 



ii,P.T. 



f8 



fcl||fc2||(fcl|| + fc2||) 



J dkidk2 PL{kl)PL{k2 

(ki -1^2)21 



fc2A:||ki+k2|2(Aq)3 



1 - 



x{sin[(fci|| -I- k2\\)Aq] — sin(fci||Ag) — sin(fc2|| Aq)}. 

(25) 

However, this choice implies that the power spectrum 
Pll(fc) only agrees with perturbation theory at the linear 
level. Indeed, it misses the contribution from the cross- 
correlation (A^ii (A\I'x)2)c and one-loop contributions 
to the variances ((Avl'||)2)^ and ((A4'_l)^)c- In princi- 
ple, this could be included by building an ansatz for the 
bivariate distribution 'P{A'i>^^, A^ ±) that does not fac- 
torize and is parameterized by the exact second and third 
cumulants (up to one-loop order). In this spirit, one may 
include perturbative results up to any finite order, at the 
price of increasingly complex models for 7'(A'I'||, A^x). 

Our approach in this paper follows a different strategy. 
We consider Eq. p9)) as a qualitative ansatz that allows 
us to build a simple model that is physically consistent 
and provides a good behavior in all regimes, as we explain 
below. Then, we treat the parameters of this ansatz as 
free parameters that we set so as to match the exact 
perturbative expansion up to the required order. In this 
paper we only go up to one-loop order (i.e.. Pi), which 

means that we choose so that Pfi^^pik) — Piioop(fc), 
that is. 



/'53(fc) =i^lloop(fc)-Piiop(fc), 



(26) 



where Pnoop is the exact one-loop power spectrum given 
by standard perturbation theory. Inverting Ea. (j24p this 
gives the effective skewness 



,eff. 



(A«) 



247r 

4 



a 



dk 



Plloop(fc)-Pliop(fc) 



2 + cos{kAq) - 3 



(Aq)4fc2 
sin(fcAg) 



kAq 



(27) 



The explicit expression p7| allows us to recover the ex- 
act one-loop power spectrum for any cosmology. To be 
meaningful, such an approach requires that we start from 
an ansatz that is not too far from the exact dynamics. 
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FIG. 3: The skewness Sg given by lowest-order perturbation 



theory, Eq. (|25|I (upper dashed line) and by the one- loop power 
spectrum matching, Ea. H27p (lower solid line). The points are 
the results from N-body simulations a.t z = 0.35 (squares) and 
z = 2 (triangles) 



This is indeed the case on a qualitative level as Eas. ([25)) 
and ((27l) show similar functional forms: in both cases 
{S^^^ (jf.^^) is a quadratic integral over P^ and S'3 " is in- 
dependent of redshift and of the amplitude of the linear 
power spectrum. On a quantitative level, the comparison 
in Fig. [3] shows that the effective skewness is greater (in 
amplitude) than the perturbative skewness by a factor of 
about two and that it has the same sign. 

Of course, if one is interested in the statistics of «;|| for 
its own sake, one should use the perturbative prediction 
([25ll rather than the effective value ([27| . In particular, 
Fig.[3]shows that the skewness is well described by lowest- 
order perturbation theory above 15/i~^Mpc, at 2 = 0.35. 
(On small scales the skewness measured in the simula- 
tions increases and becomes positive because of nonlin- 
earities and nonperturbativc effects. As expected, this 
upturn shifts to larger scales at lower redshift.) 

The effective value (|27p only makes sense as an ingre- 
dient for a model for the power spectrum. Thus, we also 
show in Fig. |4] the one-loop power spectrum (p2|) that 
would be obtained using the perturbative skewness (|25l) . 

HPT 

This yields a power spectrum P/ioop that is half-way 
between the Zcl'dovich and gravitational one-loop power 
spectra and much closer to the N-body data. There- 
fore, it already provides a more realistic starting point 
than the Zel'dovich approximation. Then, modifying 
the skewness as in Eq. ipT)) to recover the exact one- loop 
power spectrum is not a very large modification and we 
can hope that the ansatz determined by Eqs. ((T^ and 
(P7)) remains sufficiently close to the exact dynamics to 
be useful. 



7 




k [h Mpc" 



FIG. 4: Matter density power spectrum at z = 0.35, as in 
Fig. [1] We show the results from N-body simulations (data 
points), the standard one-loop prediction Pnoop, the one- loop 
prediction -Pnoop l!18|l within the Zel'dovich approximation, 
the one-loop prediction Pi{^^p' (|22|) using the perturbative 

skewness (|25p . and the full nonlinear prediction P" of Eg. psp 
using the effective skewness (|27|) . 



3. Resummed ansatz 



The one-loop power spectrum ([22|) was obtained by ex- 
panding the exponential ^ up to order P|, which only 
involves the second and third cumulants. However, one 
of the main goals of this paper is precisely not to expand 
the exponential ([9]). Indeed, it is well-known that no 
probability distribution exists such that all its cumulants 
vanish beyond some order n > 3 (i.e., the only distribu- 
tion with a finite set of nonzero cumulants is the Gaus- 
sian). In particular, approximations such as the Gram- 
Charlier or Edgeworth expansions, where we expand up 
to some finite order over cumulants, lead to functions 
that are not guaranteed to be positive and may not be 
valid probability distributions. Of course, this is not nec- 
essarily a problem, if one restricts to the range of validity 
of such approximations. However, the spirit of this pa- 
per is to avoid as far as possible such inconsistencies and 
to develop an approach that remains well-behaved in all 
regimes. 

Indeed, we can hope that by satisfying such constraints 
the convergence of our scheme will be improved and 
the matching to the highly nonlinear regime will be 
smoother. In particular, this should avoid the bad be- 
havior encountered in the standard Eulerian perturbative 
expansion, where higher orders improve the accuracy on 
very large scales but yield increasingly divergent quan- 
tities on small scales. Thus, the rise of i^Jioop at high 
k in Fig. HI around k > 0.5/i/Mpc, is an artifact due 
to the truncation at one-loop order. As for the Eule- 
rian perturbative expansion, higher orders will also give 



large contributions on these scales and the much smaller 
full nonlinear power spectrum will result from a com- 
pensation between these much larger terms (indeed, as 
explained in point ii) in Sec. IIII Al we know that the log- 
arithmic nonlinear power spectrum k^Pik) should not 
grow much beyond the nonlinear scale as long as the shell 
crossing regime has not been drastically modified from 
the Zel'dovich dynamics by trapping particles in small 
scale structures). This means that the perturbative ex- 
pansion has a very slow convergence at best (provided it 
is convergent, which is not always the case [2|), and that 
including an increasing number of higher order terms is 
not sufficient to build a useful model for our purposes. 

To resum all cumulants, it is convenient to define the 
cumulant-generating function ^p\\{y) of K|| as 



which can be expanded at y = as 



2("-l) ■ 



(28) 



(29) 



(This function is always well defined by the average i 
at least along the whole imaginary axis, even when the 
series is divergent or cumulants beyond a finite order 
do not exist.) The choice of the generating function (py 
will define our ansatz ([T9| . From the Taylor expansion 
([29| we have the behavior at the origin 



(30) 



Defining the probability distribution 'Pii(kii) through its 
cumulant-generating function ip\\ (y) shows several advan- 
tages. First, imposing the expansion ([5(1)) up to order y^ 
at the origin automatically ensures that the probability 
distribution is normalized to unity, its mean is = 1, 
its variance cr^^, and its skewness iSg". Second, choos- 
ing coefficients Sn^ (i.e., a function ip\\) that do not de- 
pend on time or the amplitude of the linear power spec- 
trum automatically ensures the scaling over Pl given by 
perturbation theory (at leading order). Indeed, as we 
have already noticed above, in perturbation theory (k|)c 

scales as P£~^- This ensures that the power spectrum 
(jl9P can be expanded over integer powers of Pl as in 
the exact perturbation theory (but of course, the exact 
value of each term will not be reproduced by our simple 
model). Third, from Ea. p8)) we can see that the aver- 
age (..)|[ that enters Eq.([T9l) has an explicit expression in 
terms of ip\\. This avoids performing an additional inte- 
gration over the probability distribution 7^||(k||), which 
is convenient for numerical purposes. In this paper we 
choose the following ansatz for V3||(j/), 



l-a 



1 - a 



(31) 
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FIG. 5; The probability distribution of the longitudinal 
relative displacement, K|| — Aa;||/Ag, a.t z = 0.35, for 
the Lagrangian separations Aq = 12.5/i~^Mpc and Aq = 
42.5/i~^Mpc. We show the linear-theory Gaussian (dotted 
line) and our model (|34p (solid lines), using for the skewness 
either Eg. pSf) (intermediate curve) or (|27|) (largest departure 
from the Gaussian). The points are the results from N-body 
simulations. 



where the scalc-dcpcndcnt parameter a{Aq) is set from 
Eq.dSni), 



a — 1 



(32) 



This gives the dependence of on the scale Aq through 

Ea. (|27l) . As seen in Fig. [3l the skewness jg nega- 

tive hence the parameter a falls in the range 



1 < a < 2. 



(33) 



The case of a zero skewness corresponds to a = 2, where 
the generating function (|3ip becomes the quadratic poly- 
nomial (y) = 1/^/2 and we recover the Gaussian case 
(i.e., the Zel'dovich approximation ([T51) ). 

The power law (pi]) is the simplest function that 
obeys the constraint (|30p while being a valid cumulant- 
generating function (this requires for instance that 
—ip\\{y) be convex). Thus, it provides a simple contin- 
uation of the series of cumulants: instead of truncating 
at third order (i.e., at the skewness) we automatically 
generate all higher-order cumulants in a manner that pro- 
vides a meaningful probability distribution V\\{k\\) in all 
regimes (i.e., 'Pii(kii) is always positive, normalized to 
unity, and with a mean = 1). From Ea.(|28p this 
"perturbative" probability distribution 'P||(k||) is given 
by the inverse Laplace transform 



^ll('^ll) = 



+ioo 



dy 



[Ky-vii(y)]/o-^ 



27ricr2, 



(34) 



Because the skewness S^, is negative, the right tail is 
sharper than the left tail, as can be checked in Fig. [5] In 
agreement with Fig. |3l Fig. [5] shows that using the per- 
turbative prediction (j25|) . with the ansatz pip , provides 
a good match to the the full probability distribution mea- 
sured on large scales in N-body simulations. As expected, 
using the larger value (P7|) for the skewness amplifies the 
deviation from the Gaussian and yields a distribution 
that is somewhat too skewed towards the left. However, 
the point of Fig.[5]is only to check that our model, defined 
by Eqs. (P7|) and ([5T|) . provides a physically consistent and 
realistic ansatz. This should be understood as a simpler 
effective system that is not expected to simultaneously 
reproduce with the same accuracy all statistical proper- 
ties of the dynamics. Thus, depending on whether one 
is interested in 7'||(k||) or in P{k) one should use Eg. ([25)1 
or (EZl). 

In both cases, the skewness is negative, which can 
be understood from the effects of gravity, as compared 
with the Zel'dovich dynamics (associated with the linear- 
theory Gaussian). Indeed, (considering for instance 
a spherical configuration) particles that move outward 
(k|| > 1) will be slowed down by gravity until they turn 
around and eventually merge into a single halo. This 
leads to an Eulerian separation that is smaller than the 
one predicted by linear theory, whence a sharper cut- 
off of the probability distribution 7'j|(K||) for large posi- 
tive In contrast, particles that move inward will be 
accelerated by gravity (because of the behavior of 
the three-dimensional gravity). This leads to a smaller 
value of the Eulerian separation than the linear predic- 
tion, which now gives a more extended tail towards low 
values of Ky, below the static value Aty, static = 1- Of 
course, this argument does not extend to k\\ < where 
shell crossing and changes of direction within virialized 
halos are expected. 

The spirit of our approach is different from resumma- 
tion schemes, where one performs partial resummations 
of higher-order terms by computing a subset of higher- 



order Feynman diagrams 
some auxiliary parameter 




or expanding over 



(such as the dimension 
of space d or the number of field components A'^). Here 
we do not explicitly compute approximations to S'„" for 
n > 4 from a subset of diagrams, as these coefficients 



are automatically generated by the functional form (j3ip 
from the low-order ones (here we only take " as input, 
but as in resummed perturbative approaches we could 
exactly include all terms up to order n and generate ap- 
proximations for higher orders by using an ansatz for t^y 
that is parameterized by the set {S'g " , .., Sn^^ }). The ad- 
vantage of our approach is that it ensures a well-behaved 
"resummed" (or rather "regularized" ) probability distri- 
bution ^[[(kii), while explicit perturbative approaches do 
not always guarantee that their resummation is physi- 
cally consistent. 

As explained in the introduction, we can hope that by 
ensuring physical consistency our method will be better 
behaved and show a faster convergence. In this approach 
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we only ensure self-consistency up to a partial degree, in 
terms of the probability distribution 7'||(A^||), but we do 
not show that physical density fields can exactly realize 
such a probability distribution, which is a more difficult 
problem. Nevertheless, this is already a significant im- 
provement over previous Lagrangian-space methods that 
expand the exponential ^ and truncate at some finite 
order [ll,!!^,^. 

Then, Eqs.(|3T|), (|32]), and ^ fully define the power 
spectrum p^ . Using Eg. ipS)) this yields 



Pll(fc) 



dAq 



(27r)3 
where we denote 



y = ikfiAqa^ 



(35) 



(36) 



We recover the Zel'dovich approximation ([T5)) for 
'/'II (y) = y ~ which corresponds to a = 2 and 

S-g" = in Eqs.dSI]) and Thus, we can sec that 

our ansatz (j35p is a continuous generalization of the 
Zel'dovich power spectrum, parameterized by the skew- 
ness S'g " . In agreement with a well-known mathematical 
result, it is only in the Gaussian case (i.e., the Zel'dovich 
approximation) that the cumulant-generating function is 
a finite order polynomial. 

As seen in Fig. 01 the power spectrum ([55]) follows the 
one-loop power (|22|) on large scales, which is identical 
to the exact one-loop power spectrum because we use 
the skewness pT]). and it improves over the Zel'dovich 
approximation (compare with the Zel'dovich one-loop 
power -Pj^oop or with Fig. [T]). At high k, P'l(fc) converges 
back to a behavior similar to the Zel'dovich power spec- 
trum, due to the fact that particles still escape to infinity 
after shell crossing. Although we will modify this behav- 
ior in the next section, to include the building of pancakes 
in a simplified form, this shows the main property that 
we looked for in this section: the power spectrum re- 
mains well behaved at high fc, with a universal behavior 
that does not change with the order of the perturbation 
theory up to which one requires consistency. 

The perturbative expansion of the power spectrum 
P5p over powers of the linear power spectrum is likely 
to be divergent for k > k"^'^^, with 



fc""^'^ = min 

Aq 



1 



(37) 



because of the finite radius of convergence of the gener- 
ating function (p\\{y) around y = 0. This gives fc™'^'^ ^ 
0.2/iMpc~^ at z ~ 0.35 (and k"^'^^ decreases with time 
as 1/(T^||). This does not necessarily apply to the actual 
power spectrum built by the gravitational dynamics and 
these values would be modified by using different gener- 
ating functions than pip. However, it explicitly shows a 
possible limitation of some perturbative approaches and 
the importance of choosing appropriate expansions or re- 
summations. This is independent of the limitation due 



to shell crossing and only due to the fast growth of high- 
order cumulants shown by our ansatz. 

The power spectrum ([35]) is closely related to La- 
grangian perturbation theory. Indeed, it can be ex- 
panded over integer powers of and it is based on a 
partial implicit resummation of higher-order cumulants 
((A^Pii )")c, which scale at leading order over as in per- 
turbation theory. For simplicity, we have only included 
the linear-theory variances ((A4'||)^) and ((A'4'x)^) and 
the third-order cumulant ((A\['||)'^)c, which is set by 
the matching with the exact one-loop contribution to 
the power spectrum. However, this approach could be 
made exact up to arbitrary order of perturbation the- 
ory by including loop contributions to these quantities 
and higher-order cumulants, such as ((A\E'||)^)c, and 
cross-correlations, such as (A^fn (A4'j_)^)c. As explained 
above, in contrast to the standard Eulerian perturbation 
theory, the resulting power spectrum would remain well 
behaved at high fc, with a Zel'dovich-like damping. 

An advantage over the Eulerian-based approaches is 
that we are not sensitive to the "sweeping effect" @,H,[23| 
associated with random advcction by long wavelengths 
of the velocity field, which only move structures by a 
random uniform shift, because we directly work with 
relative displacements. As compared with most Eule- 
rian perturbative resummation schemes, this means that 
we do not need to build a model (from a phenomeno- 
logical approach or a partial resummation of perturba- 
tive diagrams) for different-times Eulerian response func- 
tions or propagators, which are governed by the one- 
point velocity distribution 0, [13] ■ This should make such 
Lagrangian-space approaches more robust because they 
are not sensitive to such quantities that introduce addi- 
tional approximations. 



C. Adhesion-like shell-crossing continuation 

Whatever the order up to which standard perturba- 
tion theory is exactly recovered by the power spectrum 
pll(fc), following the procedure described in the previous 
section, one ingredient is still missing: the nonperturba- 
tive shell-crossing regime. In contrast to the standard 
Eulerian perturbation theory that becomes meaningless 
after shell crossing (the equations of motion themselves 
no longer make sense in the multi-streaming regime), the 
Lagrangian approach, which is based on particle trajecto- 
ries, still makes sense after shell crossing. For instance, in 
the Zel'dovich approximation the particles keep following 
straight lines after shell crossing, x(q, t) = q + '4'i(q, t). 
A similar behavior is implicitly included in the La- 
grangian ansatz psp . Shell crossing does not appear in 
this formalism and particles eventually escape to infin- 
ity as in the Zel'dovich approximation: the distributions 
P||(Aa;||) and P_l(Ax_l) widen with time. 

This implicitly provides a regularization (or comple- 
tion) of the Eulerian hydrodynamical equations of mo- 
tion. Indeed, as noticed above, the latter do not pro- 
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vide a complete description of the dynamics because they 
break down after shell crossing and one should explicitly 
state how particles behave afterwards. For instance, the 
Zel'dovich dynamics can also be written in terms of the 
Eulerian equations of motion 0, [iM HH • As compared 
with the gravitational case, one replaces the gravitational 
potential (which is coupled to the density field through 
the Poisson equation) by the velocity potential (this is 
exact at linear order). This yields an Euler equation for 
the velocity field that decouples from the density field and 
gives back the trajectories predicted by Lagrangian lin- 
ear theory. As in the gravitational case, this Euler equa- 
tion breaks down at shell crossing (when multi-streaming 
appears) and to continue the dynamics at later times 
one needs to go back to the Lagrangian interpretation in 
terms of trajectories and state that particles keep moving 
along their initial direction forever. However, this is not 
the only possible continuation. For instance, following 
the "adhesion model" introduced in [l^ , one can choose 
to make particles stick together after collisions in order 
to mimic the trapping in gravitational potential wells. In 
fact, in dimensions greater than one several continuations 
are possible. In the most convenient geometrical formula- 
tion in terms of convex hulls and Legendre transforms one 
obtains an intricate process of halo mergings and frag- 
mentations [32], nil , but it is also possible to use another 
continuation (which has no simple geometrical interpre- 
tation) where halos cannot fragment but show complex 
motions along the shock manifold [s^l ■ These two contin- 
uations and the Zel'dovich dynamics are identical before 
shell crossing and have the same perturbative expansion 
(they obey the same fluid equations in the single-stream 
regime) , they only differ after shell crossing and this only 
appears through nonperturbative terms. 

This simple example shows that we could imagine dif- 
ferent continuations of the power spectrum (j35p into the 
shell-crossing regime. Thus, Eg. ([55)) implicitly includes a 
"Zel'dovich-like" continuation but one of the main ideas 
of this paper is that this is not the best choice. Indeed, 
even though in this section we neglect halo formation, 
or more precisely the inner halo regions, we would like 
to obtain a reasonable description of the cosmic web. 
In terms of the halo model that we use in Sec. |IV] be- 
low, the one-halo term describes the inner halo shells 
and the two-halo term, which is based on the large-scale 
power spectrum that we consider here, describes the cos- 
mic web with its pancakes and filaments. As is well- 
known, such intermediate-scale structures are erased in 
the Zel'dovich approximation as particles escape to in- 
finity [2^. This has led to the "truncated Zel'dovich 
approximation" (35l . [36j . where one removes the initial 
power at high k to suppress the damping of small-scale 
structures. More generally, a Zel'dovich-like continua- 
tion, such as ([55)1 . leads to a falloff of the density power 
spectrum at high fc, due to this erasing of small-scale 
structures, and truncatingthe initial power cannot pro- 
vide more than a tail [3,[lll,[23 (e.g., for a power-law 
linear power spectrum Phik) oc fc" with —3 < n < — 1 



we have P'^{k) - A:-3-H3(n+3)/(n-Hi) j^jgj^ j^y 

This means that Zel'dovich-like continuations are not 
sufficient to describe the cosmic web. Indeed, they miss 
the formation of pancakes and filaments that give rise to 
additional power at high k. For instance, pancakes, or 
more precisely, extended sheets of zero thickness, would 
give a contribution P{k) ^ k~^ whereas infinitesimally 
thin filaments would give P{k) ^ k~^ . In the actual case 
of the gravitational dynamics with a hierarchical CDM 
power spectrum, we do not have such universal and singu- 
lar structures (pancakes and filaments can show holes of 
various sizes) and we would rather have a (multi-)fractal 
density field with a typical scaling P{k) ^ (asso- 
ciated with a correlation function ^(r) ~ j.-i-S'j |37j-[3g| 
(another modification is the finite width of the pancakes 
and filaments, which can be associated with the inner 
halo structure and is not our concern in this section). 

On a quantitative level, as recalled in the introduction, 
previous approaches that try to match perturbation the- 
ory with a halo model give too little power on intermedi- 
ate scales, where A^(fc) ^ 1. This could be traced back to 
the fact that most resummation schemes converge back 
to the linear power at high k or even show a sharper cut- 
off. From the discussion above, one explanation is that 
to match the very large scale perturbative regime, where 
shell crossing is truly negligible, to the high-fc regime as- 
sociated with inner halo regions, we need to take into ac- 
count intermediate-scale structures such as pancakes that 
give a non-negligible contribution on transition scales. 
This means that an "adhesion-like" continuation should 
be more efficient than the Zel'dovich-like continuation. 

As recalled above, in dimension greater than one, shell 
crossing is a difficult problem. In fact, in contrast to 
the one-dimensional case it is not obvious to define shell 
crossing itself as two particles coming from opposite sides 
may join circular orbits in a potential well without phys- 
ically crossing each other. In this paper, we take the 
simple criterion that was used in [l9[ to evaluate the 
impact of shell crossing on the power spectrum. If the 
Lagrangian-space to Eulerian-space mapping is potential, 
that is, x(q, t) = i9$/9q, as in the Zel'dovich and ad- 
hesion models, shell crossing is associated with the loss 
of convexity of the potential $(q, t). (In the Zel'dovich 
dynamics, is of the form |qp/2 -I- 0^, where is 
stochastic and given by linear theory, and loses its 
initial convexity as (f>L grows with time, while in the ad- 
hesion model we take its convex hull, conv($^), which 
prevents shell crossing but gives rise to shocks |3ll - [33| ). 
Then, convexity of <I>(q, t) implies that Aa;|[ > (i.e., the 
longitudinal Eulerian separation docs not change sign, 
for any pair of particles) [iO, Hlj . Then we can choose 
Aa;|| < as a criterion of shell crossing [l^. This is ob- 
vious in one dimension and provides the adequate gen- 
eralization to higher dimensions for potential mappings. 
In the case of the gravitational dynamics, the mapping is 
only potential up to second order of perturbation theory 
[43 | , but we can still expect that it gives a useful criterion 
for the formation of the first large-scale structures. 
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FIG. 6: Skewness 5*3 given by low-order perturbation theory 



(|25|l (upper dashed line) and by the effective value (|27|l (lower 
dashed line). The solid lines that arise from these two large- 
scale asymptotes are the predictions at redshifts z — 0.35, 1, 2 
and 3 (from right to left) obtained from the nonperturbative 
adhesion-like probability distribution p8|l , determined by the 
corresponding perturbative part Py . The points are the re- 
sults from N-body simulations at 2 = 0.35 (squares), z = 1 
(diamonds), z = 2 (triangles), and z = 3 (circles). 



Therefore, we modify the "perturbative" longitudinal 
probability distribution ([M)) by setting Aa;|| = to all 
pairs that had Aa;|| < 0. More precisely, we define the 
"adhesion-like" longitudinal probability distribution as 

7'|f-(K||) = aie(AC|| >0)7'||(At||) + ao(5i5(K||), (38) 

where 6(k|| > 0) is the Heaviside function (1 for K|| > 
and for K|| < 0). The parameters ao and oi are set 
by the constraints (1) = 1 (i.e., the probability distri- 

= 1. From the 



bution is normalized to unity) and 
expression wc obtain 



ai = {1 + Ai) ^ and oo 
where we introduced 

r.O++ioo 



1^0, 



(39) 



A-n — 



dy 



n+l 



0+- 



27ricr2 



y 



-Mv)/-.,^ (40) 



where the contour over y runs to the right of the pole 
at the origin and to the left of the branch cut at = 
a — 1. The coefficients An are nonperturbative and scale 



as e 



(fli - 1) - ao 



, which gives 



) 



(41) 



Therefore, the "perturbative" distribution V\\ obtained 
in Eq.dSH) in Sec. IIII B 31 and its "adhesion- like" modi- 
fication are identical to all orders of perturbation 
theory. 



We show in Fig. [6] the skewness " defined by this 
new probability distribution psp . using either the per- 
turbative skewness (|25p or the effective value (P7)) for 
the perturbative part Vw . On large scales we recover the 
skewness associated with the regular distribution (|34p . 
with a very fast convergence because the deviation de- 
cays as ~ e"''^'-' On small scales, the skewness in- 
creases and becomes positive in a fashion similar to the 
behavior measured in the simulations. As for the neg- 
ative sign in the perturbative regime that we explained 
in Sec. IIII Bl this can be understood from the dynam- 
ics. Indeed, because of the "sticking" of particle pairs at 
Axy = 0, in the nonlinear regime (i.e., on small scales), 
which becomes sensitive to this modification, the 1ow-K|| 
tail becomes very sharp (it is zero for «;|| < 0) whereas 
the high-K|| tail still extends to infinity, albeit with an 
exponential-like decay. Therefore, the global shape of 
the probability distribution now shows a broader right 
tail, in contrast with the perturbative regime shown in 
Fig. [51 which now leads to a positive skewness. In the ac- 
tual gravitational case, there is no such exact "sticking" 
to Axy — but there is a trapping within small virial- 
ized halos. Thus, collapsed pairs remain bound with a 
separation set by the typical size a;haio of virialized ha- 
los. This plays the role of the left-tail cutoff, which is no 
longer sharp at Axy = but decays on a scale of order 
s^haio, whereas the right tail still extends to infinity and 
is not strongly affected by smaller-scale virialization (it 
corresponds to rare voids). Of course, we cannot expect 
the simple model psp to provide an accurate prediction 
for S'g " , even when we use the correct perturbative limit 
on large scales. However, we can check in Fig. [5] 
that it already provides a good qualitative description. 
In particular, it captures the dependence on redshift of 
the upturn of S'g " , due to these nonperturbative effects 
that occur after shell crossing. 

Then, we modify the perturbative power spectrum 
([55)1 by replacing the perturbative probability distribu- 
tion ([M)) by its adhesion- like extension psp . Substituting 
into Ea. (jl9p gives the "cosmic web" power spectrum 



-few. (k) 



0++ioo 
/0+-ioo 

X e 2' 



dAq 1 [ -v\\(-iHi.Aqal)/o 

(27r)3 1 + Ai^^ 

dy -v\\{y)/'yl 
27ri 



1 



1 



y y + ^k^Aqal 



(42) 



where the contour over y again crosses the real axis in- 
between 0<y<Q!— 1. As in [l^, this is a "sticky 
model" that can be seen as a very simplified version of 
the full 3D adhesion model. Since we do not modify the 
transverse motion, the "adhesion" only takes place along 
the longitudinal direction, which also serves as the signal 
of shell crossing. Therefore, we only include planar struc- 
tures (thin pancakes). To describe filaments we should 
also include some sticking along one transverse direction, 
but this would require additional ingredients and free pa- 



12 




k [h Mpc" 



FIG. 7: Matter density power spectrum at z = 0.35, as in 
Figs. [T]and|4l We show the results from N-body simulations 
(data points), the nonlinear Zel'dovich power spectrum 
(|16|l . the nonlinear ansatz P" (|35}, and the "adhesion-like" 
continuation Pc.w. ()42p for the "cosmic web". 



ramctcrs (e.g., to set the relative mass between filaments 
and pancakes). Hence for simplicity we only take into 
account pancakes as in Ea. (|42|) . which fits in a natural 
fashion in our framework where we have already kept 
track of the longitudinal separation. 

As explained above, the shell crossing correction is 
nonperturbative and the power spectrum (j42p is identical 
to the power spectrum (|35|) to all orders over Pl. 

We compare in Fig. [7] the three Lagrangian-space 
power spectra that we have obtained, the usual 
Zel'dovich approximation P^, our nonlinear ansatz P\ 
and its "adhesion-like" continuation Pc.w.- On large 
scales the nonperturbative correction is negligible and 
we recover the perturbative power spectrum P", which 
is somewhat larger than the Zel'dovich power spectrum 
because of the nonzero skewness " that allows us to en- 
sure consistency with standard perturbation theory up 
to one-loop order. On small scales, the nonperturba- 
tive correction becomes dominant and gives some extra 
power, associated with the formation of pancakes, with 
a high-fc tail ^ k^^ that decreases more slowly than the 
Zel'dovich-like tails (that are steeper than fc~^). The 
nonlinear power spectrum Pc.w. is the Lagrangian ansatz 
that we use in this paper to describe the "cosmic web" . 



IV. COMBINING THE "COSMIC WEB POWER 
SPECTRUM" WITH THE HALO MODEL 

A. The halo model from a Lagrangian point of view 

Because our goal is to build an unified model for the 
power spectrum that applies from linear to highly nonlin- 



ear scales, we must combine the perturbative approach 
described in the previous section (which is systematic 
and accurate but only applies to large scales) with phe- 
nomenological models (that are not systematic and only 
show an accuracy of 10% but can be applied to small 
scales) . Following [ll|, [l5[ , we consider the halo model 
from a Lagrangian point of view instead of the usual Eu- 
lerian framework [43| . This provides a convenient frame- 
work to include our Lagrangian perturbative approach 
within the halo model (the latter being mostly used to 
describe small highly nonlinear scales). In particular, the 
transition from the perturbative large scales, driven by 
bulk flows, to the inner halo regions, driven by virial- 
ization, takes place in a gradual fashion while satisfying 
matter conservation (i.e., without double counting). 

Following [ll|, we split the average in Eq.ijB]) over two 
terms, Pih and P2H, associated with pairs {qi,q2} that 
belong either to a single halo or to two different halos. 



P(fc) = PlH(fc)+P2H(fc), 



(43) 



with 



Pm{k) - / ^ ^iH(Ag) (e"^-^- - e"-^'>)iH (44) 
and 

^2H(fc) = / ^ ^2h(A<z) {e'^-^- - e'''-^'i)2H. (45) 

Here the averages (...)ih and (...)2H are the conditional 
averages, knowing that the pair of length Aq belongs to 
a single halo or to two halos, while Pih and P2H are the 
associated probabilities. The approximation of spherical 
halos in Lagrangian space gives for the probability to 
belong to a single halo [Tl| . 



Pih(A(7) 



(2gA/- Ag)2(4(7M + Ag) 



I6q 



(46) 

where qm is the Lagrangian radius associated with the 
mass M, M = Anjiqlj/S, and f{iy) is the scaling function 
that defines the halo mass function. 

Here a{M) is the rms linear density contrast at scale M 
and 6l = J^~^(200) is the linear density contrast asso- 
ciated with the nonlinear density threshold that defines 
collapsed halos, which we choose equal to 200. In numer- 
ical computations we use the fit to the halo mass func- 
tion from 4^ , which has been shown to match numerical 



simulations while obeying the asymptotic large-mass tail 
f{iy) ~ e-"^/'^ m. The probability that the pair belongs 
to two different halos reads as 



P2H(A<z) = l-PiH(Ag). 



(48) 



This automatically avoids any double-counting in the de- 
composition (j43|) and it ensures that we count all matter 
once. 
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B. One-halo term 



Two-halo term 



The decomposition (|43|) also corresponds to the one- 
halo and two- halo terms of the usual halo model [i^ . In 
particular, Eqs.gl]) and ge]) yield fu\.ll^ 



Av M 
— j{v)—— 
P(27i-)- 



PMk) = / —I{v)^f^ [uM{k) - W{kqM) 

(49) 

Here UM{k) is the normalized Fourier transform of the 
halo radial profile, 



(50) 

where pm{x) is the halo density profile for a halo of mass 
M, and W{kq) is the normalized Fourier transform of the 
top hat of radius q, 



sm{kq) — kqcos{kq) 



(51) 



To derive Eq. we used the approximation of fully viri- 
alized halos: the two particles qi and q2 have lost all 
memory of their initial locations and are independently 
located at random within the halo. As in [ll|, in nu- 
merical computations we use the usual NFW halo profile 



The counter term W in Eq. p9)) arises from the counter 
term e'^ '^'i of Eq.(|33]) (this subtracts the contribution 
from the mean density). The computation in [ij would 
rather give a factor {u\j — W^), but we prefer to use the 
factor (um ~ W)^ that readily extends to higher order 
multi-spectra as seen in [l^. Moreover, this ensures that 
the one-halo contribution to the matter power spectrum 
decays as fc^ at low k, as implied by the conservation 
of matter and momentum for small-scale redistributions 
of matter (whereas the factor {u\j — W"^) only en- 
sures a k^ tail, which is consistent with the conservation 
of matter but not of momentum). We show in Fig. [20] 
in App. 1^ the impact of this \ow-k tail of the one-halo 
term on the power spectrum. We find that using a fac- 
tor {u\j — W'^), which gives a slower falloff at low fc, 
can overestimate the power spectrum by about 5% on 
transition scales. Indeed, at higher k, in the highly non- 
linear regime, the counter term W is negligible, whereas 
at lower k the onc-halo term itself is negligible. This 
shows that on transition scales, where A^(fc) ^ 1, the 
power spectrum is sensitive to details of the halo model 
if we require an accuracy of a few percents. Fortunately, 
this only appears on a limited range of scales and this 
has no impact on the perturbative scales. (However, it 
is important to ensure that the one-halo term decays at 
least as k^ at low k rather than converging to a constant 
as in the usual prescription without any counter term.) 



At a perturbative level Fm and Pm are identically 
zero while F2H is unity, because of the exponential de- 
cay of the halo mass function, of the form e~^^'^ i^) ^ in 
the rare event limit. Then, as noticed in [Tl| . the power 
spectrum given by perturbation theory is included in the 
two-halo contribution (|45|) . More precisely, the expan- 
sion over powers of Pl of P2H must recover the standard 
perturbative expansion. In [llj, we used the simple ap- 
proximation P2h(^) — P2H(l/fc)Ppert(fc), whcrC Ppcrt(fc) 

is the power spectrum given by the Eulerian "steepest- 
descent" resummation scheme developed in 0, . Any 
other Eulerian or Lagrangian resummation scheme could 
be used, provided it is well behaved at high k where it 
becomes subdominant with respect to the one-halo term 
(this excludes the standard perturbation theory, which 
grows too fast at high k, unless one adds an extra high-fc 
cutoff). 

However, the approaches investigated in [ll| did not 
manage to provide a fully satisfactory matching to the 
highly nonlinear regime because they predicted too lit- 
tle power on intermediate scales (where A^(A:) ~ 1 — 10). 
This can be traced to the fact that they usually go back to 
the linear power at most at high fc, which leads to insuf- 
ficient power on transition scales where P2H and Pih are 
of the same order. Another problem was that Lagrangian 
perturbative schemes, which would be more convenient 
to embed within Eq. ps)) and to extend to redshift space, 
made this lack of power even worse because they usually 
display a strong cutoff at high k. In this paper, our goal 
is to improve over [T]| by implementing a Lagrangian 
perturbative scheme that is free of this problem and pro- 
vides reasonably accurate predictions up to the transition 
scales. Thus, we use the "cosmic web" power spectrum 
developed in Sect. [Hi] as a basis of our two-halo term, 
which we write as 



P2H(fc) 



(2^ ^^2h(A<z) (e 



_i 1.2 
X e 2 



0++ioo 



0+-ioo 



dy -vii(y)/<^2 
27ri ^ 



1 



1 



y y + ikp.Aqal 



(52) 



We recognize the cosmic web power spectrum given 
by Ea. (|^^ to which we added the factors P2H and 
^gik-Ax^vir ^ As explained in Eqs.(l43l)-(l45]), the factor P2H, 
given by Eqs. ipS)) and (|46p . ensures that mass is con- 
served: there is no double counting of particle pairs. The 
factor (e'*''^'')^^ is associated with small-scale motions. 

Indeed, as seen from the derivation of Eq. p^ in 
Sec, mil the cosmic web power spectrum Pc.w. arises from 
the statistical average of eJ^'^^ due to large scale mo- 
tions, associated with the bulk fiows described by pertur- 
bation theory that we regularize by an adhesion-like con- 
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tinuation. As noticed in Scc. lIII Cl this simple procedure 
takes into account the formation of idealized infinitesi- 
mally thin pancakes. In the actual gravitational process, 
pancakes and filaments have a finite width, as particles 
keep moving for some time after shell crossing instead 
of instantaneously sticking together and have finite-size 
turn-around radii and virialized orbits. More generally, 
we can split the motion of particles in two components, a 
first one associated with large-scale bulk fiows, which was 
considered in Sec. Illll and corresponds to the "skeleton" 
of the large-scale structures, and a second one associ- 
ated with small scale virialized motions, which gives some 
thickness to this skeleton. In particular, for hierarchical 
linear power spectra with a high-fc tail that does not de- 
crease faster than all particles are expected to belong 
to collapsed objects (this may not be the case for linear 
power spectra with less power at high k). Then, making 
the approximation that these small-scale and large-scale 
motions are decorrelated, we write 



ik-Ax 



ik-Ax\bulk / ik-Ax\vir 



(53) 



The first part was the focus of Sec. IIIII and corresponds 
to the cosmic web power spectrum Pc.w.- Assuming as 
in the derivation of the one-halo term (|49p full virializa- 
tion [^], that is, that the two particles qi and q2 are 
independently located at random in their host halo, we 
write 



ik- Ax\ vir 



^glk-ZXx^ 



'Ah 



(54) 



In contrast with the one-halo case where the two 
particles belong to the same halo, which gave rise to the 
factor uj^, here the two particles belong to two different 
halos whence the integration over two halo mass func- 
tions, each one with its factor um- In Eq. (|54P we only 
integrate over halos of radius smaller than Aq/2, to take 
into account in an approximate fashion that if a pair of 
separation Aq belongs to two different halos these halos 
are unlikely to have a radius much greater than Aq/2. 
Thus, massive halos only contribute to the one-halo term 
(Uni) and to the two-halo term ([5^ at large distances. 

We show in Fig.l^Olin App.|X]the impact of this "virial 
damping" factor on the matter power spectrum, by plot- 
ting the deviation that would be obtained by neglecting 
this term. Because of the upper bound on halo mass in 
Eq. (|54p , this factor does not induce a significant damping 
of the power spectrum ((52)) on large scales. Indeed, in the 
weakly nonlinear regime where the two-halo term is dom- 
inant, power at wavenumber k typically comes from pairs 
of initial separation Aq ~ 1/fc, which selects halo radii 
qM that are smaller than 1/fc for which the factor UM{k) 
is close to unity. On small highly nonlinear scales it yields 
a greater damping as 1/fc becomes smaller than the typi- 
cal size of the halos. However, on these scales the matter 



power spectrum given by our approach is not accurate 
to better than 10% because of the one-halo term itself, 
which involves the halo profiles and their concentration 
parameters that are not modeled to a very high accuracy. 
Nevertheless, we include this factor in our computations 
because it naturally arises in our framework. 



D. Nonperturbative effects on large scales 



Thus, in our approach we include nonperturbative and 
shell crossing effects on the large scale power spectrum 
through the factors F2H and (e''''^'')^^, associated with 
halo formation and virialized motions, and through the 
adhesion-like continuation described in Sec. IIII CI (in ad- 
dition to the one-halo term itself). This is different from 
recent Eulerian-space works 4^, 4^ that propose to in- 
clude the impact of such nonperturbative effects on large 
scales through additional terms to the hydrodynamical 
equations of motion of the fluid approximation, that may 
be obtained by a combination of coarse-graining (that 
sets the form of these terms) and phenomenology (they 
include coefficients that are measured in simulations). In- 
deed, in our approach the hydrodynamical equations of 
motion enter through the perturbative expansion that 
they imply for the power spectrum, which is not affected 
by these nonperturbative effects. Then, the impact of a 
small-scale velocity dispersion, due to virialized motions. 



is described by the factor (e''' '^^)^', in Eq.dSS]). This 
would play the role of some pressure-like terms included 
in [4^, The trapping within potential wells, that 

is described in a simplified manner by the nonpertur- 
bative correction associated with the adhesion-like 
continuation, may also be described within such methods 
through a pressure or viscosity term, as in the original 
adhesion model [Tsl where the pressureless Euler equa- 
tion is replaced by the Burgers equation. 

Another approach to take into account such effects 
would be to go back to the Vlasov equation of motion (sol - 
[52| . In principle, this could provide systematic schemes 
but the methods that have been proposed so far lead to 
heavier computations and no fast and accurate method 
has been presented yet. 

As noticed in the introduction, the advantage of 
the Lagrangian-space framework over the Eulerian-space 
framework is that particle trajectories can describe both 
the single-stream and multi-stream regimes, which allows 
us to include these nonperturbative effects in the simple 
fashion described in the previous sections. However, this 
approach is not fully systematic and it involves some phe- 
nomenological insight. 



15 



V. NUMERICAL RESULTS 

A. WMAP5 cosmology 

We first compare our model with numerical simulations 
for a WMAP5 cosmology, that is, with the set of cosmo- 
logical parameters h = 0.701, 17,n = 0.279, as = 0.8159, 
Us = 0.96, and Wdc = —1- While we use the N-body 
data of [iJI on large scales (fc < 0.3h Mpc~^), we use 
the measurements of [llj on smaller scales. This is be- 
cause the former one is tuned to give an accurate power 
spectrum on BAO scales and has a larger total volume 
(60 X (2048 Mpc)3 - 515 h'^ Gpc^), while the latter 
gives more reliable results over a wide range in k by care- 
fully combining the measurements from some simulations 
with different resolutions (the highest-resolution simula- 
tion employs 2048^ particles in a (512 ft,^^ Mpc)'^ cube). 
Similarly, we use the simulation results from these papers 
for the two-point correlation function: data points from 
[11] are shown on BAO scales (from 70 to 140 Mpc), 
while those from [lll | are depicted on smaller separations. 

1. Power spectrum 

We show in Fig. [S] the final power spectrum obtained 
by our model, combining the cosmic web power spectrum 
obtained in Sec. IIIII with the halo model. We clearly see 
how the two-halo and one-halo contributions arc domi- 
nant on large and small scales respectively. The transi- 
tion takes place around k ~ lft,Mpc~^ at low redshift, 
in a gradual manner. In contrast with previous studies 
[ill [isl l , we no longer underestimate the power spectrum 
on these transition scales and we obtain a good agreement 
up to the highly nonlinear regime. This is because our 
two-halo contribution is no longer based on a perturba- 
tivc power spectrum that goes back to the linear power 
at high k but on the "cosmic web" power spectrum of 
Sec. mil which shows more power on nonlinear scales be- 
cause it takes into account intermediate-scale structures 
such as pancakes. This gives a more realistic basis that 
appears indeed to provide a better match to the N-body 
results. 

The upper row in Fig. |5] clearly shows how the slope 
of the power spectrum evolves with redshift. For CDM 
power spectra, the transition scale shifts to higher k at 
higher redshift, which leads to a "redder" power spec- 
trum (the local slope n of the linear power spectrum de- 
creases and typically shifts from —1.5 to —2.5). In turn, 
this yields a two-halo contribution that becomes larger 
with respect to the one-halo contribution on highly non- 
linear scales. In retrospect, this explains why the un- 
derestimation of the power spectrum on transition scales 
was more severe at higher redshift in [llj, [l^ . Although 
using the "cosmic web" power spectrum of Sect. |lll]is a 
significant improvement over these previous works, our 
two-halo contribution is unlikely to be very accurate at 
high k and this is probably one of the reasons for the 



discrepancy found in the rightmost panels of Fig. [8] at 
k > 3/iMpc^^. On these small scales, we may overesti- 
mate the two-halo contribution. 

Another source of inaccuracy in this regime is the un- 
certainty of the halo-profile parameters themselves, in 
particular the mass-concentration relation c{M). We 
show in Fig. [5] the results obtained by using for c(M) 
the fits to N-body simulations from [5J| and f54l|, as well 
as the formula 

/ M \""-^ 

for 0.35 < z < 3, where halos are defined by a den- 
sity contrast of 200 with respect to the mean density of 
the universe. This gives a power spectrum that is in- 
between those obtained with the fits from [s^ and [HI] 
and gives a slightly better match to our simulations at 
high k. In principle, the relation c{M) used for the power 
spectrum is expected to be slightly different from the 
one measured in simulations because of its finite scatter, 
as inclusion in the matter power spectrum is associated 
with an implicit weight that may differ from the aver- 
aging procedure used to measure the relation c{M) in 
the simulations. In any case, we can see in Fig. [8] that 
the discrepancy between our model and the power spec- 
trum measured in the N-body simulation is of the same 
order as the difference between the predictions that use 
either [5^ or for c(M). This means that the power 
spectrum at these high wavenumbers still shows an un- 
certainty on the order of 10%, whether it is computed 
through the halo model or directly from the N-body sim- 
ulations. In particular, shot noise certainly explains part 
of the rise of the power spectrum measured in the N- 
body simulations at z = 3 for k > 20/iMpc~"'^ above our 
predictions (lower right panel in Fig. ^ . This is clearly 
seen from the comparison with the panels in the lower 
row obtained at lower redshift, which probe deeper into 
the nonlinear regime before being affected by shot noise, 
where we can see that the logarithmic power spectrum 
follows the shape predicted from the halo model. In par- 
ticular, at z — 0.35 we clearly see the slowing down of 
the growth of A^(fc) — Ank^P(k) in the highly nonlinear 
regime and we would expect a similar behavior at z = 3 
(actually we would expect a slightly faster slowdown be- 
cause the local slope n of the linear power spectrum is 
redder). Therefore, in this regime it seems that semi- 
analytical approaches like ours, based on the halo model, 
are competitive with direct N-body simulations. (At very 
high k the semi-analytic approaches are expected to re- 
main reasonable, because they are based on a physically 
reasonable ansatz and/or assumption, whereas the direct 
results from simulations suffer from shot noise). 

We show in Figs. [9] and [10] the relative deviation be- 
tween the power spectrum predicted by our model and 
the N-body measurements. We obtain an accuracy of 
about 2% up to A; '--^ 0.3/iMpc~^, and 5% up to fc '--^ 
3/iMpc~-'^, for z > 0.35. The accuracy degrades rapidly 
in the highly nonlinear regime, because of the uncertainty 
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FIG. 8: Upper row: matter density power spectrum P{k) at z = 0.35, 1, and 3, from left to right, multiplied by a factor [kh)^'^ 
as in Fig. [7] Lower row: logarithmic matter power spectrum, A'^(fc) = 4nk^ P(k). We show the linear power spectrum, the two- 
halo and one-halo contributions P2H (|52|) and Pm (|49|l . the full nonlinear power spectrum (|43p . and the results from numerical 
simulations. For the full nonlinear power spectrum (|43p . in addition to the result obtained with the mass-concentration relation 
(I55|l (solid line) we also show the results obtained using the mass-concentration relations from [H^l and [sj (dashed lines). 



of the halo model and of the N-body simulations them- 
selves (in particular because of shot noise at very high 
k). Fortunately, as shown in Fig. [HI the uncertainties 
of the halo model (i.e., the parameters of halo profiles) 
do not contaminate the predictions for the power spec- 
trum on large scales, k < l/iMpc~^. Therefore, these 
large scales remain a robust probe of cosmology. This 
is also one interest of such analytical approaches that is 
complementary to numerical simulations: they allow us 
to estimate the impact of different processes on the final 
power spectrum and to estimate the range of wavcnum- 
bers that is not affected by small-scale uncertainties and 
can be safely used to constrain cosmology up to a good 
accuracy. 



2. Two-point correlation function 

Wc show in Fig. [TT] our results for the matter density 
two-point correlation ^(x), given by 

ax) = in r dfc eP{k) 5^^. (56) 

As in previous works [ll|, [l^ [s^ [s^, we can see that 
using a well-behaved perturbative contribution that in- 
cludes one-loop contributions provides a good accuracy 
on baryon acoustic oscillation (BAO) scales. In par- 
ticular, we reproduce the well-known damping of the 
baryonic peak at ~ 105/i^^Mpc as compared with lin- 
ear theory. We can see in the lower panel the transi- 
tion between the two-halo and one-halo contributions, at 
X ~ l/i~^Mpc. As for the power spectrum, we obtain 
a significant improvement over previous studies (ill . [l5j 
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FIG. 9: Relative deviation between the density power spec- 
trum predicted by our model and the result from numerical 
simulations, at z — 0.35 (black squares), 1 (red circles), and 3 
(blue triangles) . We show the results (which almost coincide) 
obtained using the mass-concentration relation (|55p and the 
ones from ^31 and [5411. 



Q. 

2~ 




0.1 1 10 

k[h Mpc"^] 



cause the local slope of the linear power spectrum on 
transition scales becomes redder while our "cosmic web" 
power spectrum shows a universal fc~^ tail. The good 
match with simulations on transition scales suggests that 
using such a two-halo contribution that is greater than 
the linear one on nonlinear scales is an important ingre- 
dient to recover the full nonlinear power spectrum. 

On small scales, the discrepancy between our model 
and the numerical simulations is again on the order of the 
scatter due to the uncertainty of the halo-model mass- 
concentration relation. Fortunately, this only affects very 
small scales, where the physics of baryons should also be 
included (for instance, through the halo-model parame- 
ters) [st} . 

We show in Figs. [T^ and [T51 the relative deviation be- 
tween the two-point correlation predicted by our model 
and the one measured in the N-body simulations. We 
obtain an accuracy of about 2% on the BAO scales and 
of about 5% down to O.S/i^^Mpc. A gain, the accuracy 
degrades on smaller scales because of the uncertainties 
of the halo-model parameters and of the finite resolu- 
tion of the N-body simulations themselves. On very 
large scales the statistical error bars of the simulations 
grow because of the finite-box size and become much 
larger than the inaccuracy of the analytical model. In 
fact, it seems that we predict more power than is mea- 
sured in the simulations by about 2% on large scales. 
Part of this offset may be due to a lack of large-scale 
power in the N-body simulations because of their finite 
size. Indeed, we obtain a better agreement on interme- 
diate scales, 20 < a; < 60/i~^Mpc. Moreover, the two- 
point correlation function changes sign and vanishes at 
xq ^ 130/i^^Mpc, so that small absolute deviations are 
amplified after we take the ratio to 'Csim(a;) and the rel- 
ative deviation becomes infinite at xq because of the fi- 
nite accuracy of models and simulations. Thus, as for the 
power spectrum, semi-analytical approaches like ours ap- 
pear competitive with direct N-body simulations. 



B. Other cosmologies 



FIG. 10: Relative deviation between the density power spec- 
trum predicted by our model and the result from numerical 
simulations, at z — 0.35 (squares), 1 (circles), and 3 (trian- 
gles). The solid lines are estimates of the effect of the shot 
noise in numerical simulations (smaller impact at lower red- 
shift). 



To further test our model, we also compare our predic- 
tions with numerical simulations for two other cosmolo- 
gies. 

1. WMAP3 cosmology 



as we no longer underestimate the two-point correlation 
on these transition scales. In particular, we recover the 
shape of the two-point correlation from linear to highly 
nonlinear scales. Again, the two-halo contribution be- 
comes significantly greater than the linear correlation 
(or the linear power spectrum in Fig. [5]) on small scales 
because of the adhesion-like continuation explained in 
Sect, nil CI This is more important at higher rcdshift be- 



We first use an older set of WMAP3 simulations, with 
h = 0.734, fJ„i = 0.234, CT8 = 0.76, = 0.961, and 
ifdo = ^1, performed in [58|. These simulations have a 
smaller box size (1000 Mpc) and a lower resolution 
(512'^ particles), hence we only perform the comparison 
on large scales. On the other hand, they allow us to go 
down to redshift 2 = 0. 

We show our results for the power spectrum in Figs. 1141 
and [ini These results are similar to those obtained in 
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FIG. 11: Upper row: matter density correlation function ^(a::) at z = 0.35, 1, and 3, from left to right, multiplied by a factor 
10^, on BAO scales. Lower row: correlation function ^(x) on smaller scales. We show the linear correlation ^l, the two-halo 



and one-halo contributions ^2h and ^ih from Eqs. (|52|) and (|49|) . and the full nonlinear correlation ^ from Eq. (|43p (using the 
mass-concentration relation (|55|l (solid line) and the fits from [S^l and [sj (dashed lines)). The points are the results from 
numerical simulations. 



Figs. H] and |9] for the WMAP5 cosmology, with larger 
error bars in Fig. [15] because of the smaller box size of 
the simulations. 



2. Quintessence cosmology 

We finally consider a less standard cosmology, where 
the dark energy does not behave as a cosmological con- 
stant. Thus, we choose a dark energy equation of state 
parameter u>do ~ —1.281, and h = 0.7737, fim = 0.23638, 
(78 = 0.7692, and Us = 1.0177. We use the measurements 
from N-body simulations done in [l^ for this cosmologi- 
cal model. These simulations are done in 2048 ft." ^ Mpc 
boxes with 1024"^ particles, and eight statistically inde- 
pendent realizations are available. 

We show our results for the matter power spectrum 
in Figs, [in] and 1171 and for the correlation function in 



Figs. [18] and [19] We again obtain similar behaviors 
to those found for the LCDM WMAP5 cosmology in 
Sect.[V3 

This confirms that our method provides efficient pre- 
dictions that can be used for a variety of cosmologies. 

VI. CONCLUSION 

Wc have developed in this paper a new approach to a 
systematic modeling of the cosmological density field: in- 
stead of looking for explicit partial resummation schemes 
we embed the standard perturbation theory within a re- 
alistic nonlinear ansatz. This automatically ensures a 
reasonable behavior on small scales, while being consis- 
tent with perturbative approaches up to the required or- 
der. Then, this allows us to obtain a reasonable matching 
with the highly nonlinear regime. 
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FIG. 12: Relative deviation between the two-point correlation 
^(x) predicted by our model and the result from numerical 
simulations, at 2 = 0.35 (squares), 1 (circles), and 3 (trian- 
gles). We show the results (which almost coincide) obtained 
using the mass-concentration relation (|55p and the ones from 
[5^ 1 and [131. The error bars are the statistical error bars of 
the N-body simulations. 
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FIG. 13: Relative deviation between the two-point correla- 
tion ^(x) predicted by our model and the result from nu- 
merical simulations, at z = 0.35 (squares), 1 (circles), and 3 
(triangles). The error bars are the statistical error bars of the 
N-body simulations. The solid lines are estimates of the effect 
of the shot noise in numerical simulations (smaller impact at 
lower redshift). 



We have shown how this can be achieved within 
a Lagrangian-space framework. Because the lowest- 
order approximation that app ears in this context is the 
Zel'dovich approximation [23|, we first discussed the lack 
of power on weakly nonlinear scales associated with this 
Gaussian approximation. Then, we explained how this 



can be partially cured by including higher-order terms 
(i.e., the skewness) while keeping a well-behaved ansatz. 
Because the Gaussian is the only probability distribution 
with a finite set of nonzero cumulants, this requires in- 
cluding nonzero cumulants at all orders and we describe 
a simple ansatz that provides a well-behaved probabil- 
ity distribution for the longitudinal displacement field. 
This provides a perturbativc power spectrum that has 
the same qualitative properties as the true gravitational 
dynamics power spectrum. Moreover, it can be made 
consistent with the exact perturbativc expansion up to 
the required order, while keeping all these qualitative 
properties. Here we only go up to order P|, but higher 
orders could be taken into account by explicitly choosing 
the kurtosis and higher-order cumulants of the displace- 
ment field (at the price of more complex ansatz for the 
cumulant generating function). 

Next, we argued that for a good description of inter- 
mediate scales it is important to include some effects 
associated with the shell crossing regime. In particu- 
lar, we point out that the small-scale behavior implied 
by the Zel'dovich power spectrum and most previous 
Lagrangian-space schemes, associated with the escape of 
particles to infinity, is not adequate and a more realistic 
picture is provided by the adhesion model. We present 
a simplified sticking model that captures some features 
of this trapping of particles within potential wells and 
describes the first stages of pancake formation. In par- 
ticular, it captures the (nonperturbative) increase and 
change of sign of the skewness of the displacement field 
on small scales. This gives a "cosmic web" power spec- 
trum that shows a broader high-A: tail than the original 
Zel'dovich power spectrum, with a universal k~'^ decay 
instead of the steeper than falloff. The final expres- 
sion for this "cosmic web" power spectrum is summarized 
as [see Eq. 
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with the quantities ctj 
spcctively given by Eqs. (jH)) . and PO)) . For the 

cumulant-generating function (^||, we adopt the simple 
ansatz of Eq. ((3T|) . which is parameterized by the effec- 
tive skewness ^7} through the scale-dependent parame- 
ter a(Ag) [Eq. ([5^]. 

This provides a well-behaved cosmic web power spec- 
trum that can be combined with a halo model to describe 
all regimes of the cosmological density field, from large 
linear to small highly nonlinear scales. The explicit ex- 
pression for the power spectrum, given as the sum of 
the two contributions Pmik) and P2H(fc), is described in 
Sec. IV] [see Eqs. ^ and i^]. 

A comparison with numerical simulations shows that 
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FIG. 15: Relative deviation between the density power spec- 
trum predicted by our model and the result from numerical 
simulations, at 2: = (black squares), 1 (red circles), and 3 
(blue triangles), as in Fig. [9] but for a different cosmology 
(WMAP3). We show the results (which almost coincide) ob- 
tained using the mass-concentration relation (|55p and the ones 
from [sl and [HI. 



our new proposition extends the range of validity of the- 
oretical predictions over previous models. We obtain an 
accuracy better than 2% for k < 0.3/iMpc~^ and bet- 
ter than 5% for k < ShMpc'^, at z > 0.35. This also 
yields an accuracy of 2% on BAO scales for the two-point 
correlation function and of 5% down to 0.5/i~"'^Mpc. We 
checked that a similar accuracy is obtained for an al- 
ternative less standard cosmology, with a dark energy 
equation-of-state parameter Wdc — —1.28. 



Our approach should be generalized in two directions 
to provide a more direct comparison with observations. 
First, we should take into account redshift-space distor- 
tions, which should be possible within our Lagrangian- 
space framework. Second, we could consider the statis- 
tics of biased tracers, such as galaxies of different mass 
or luminosity. We leave these generalizations to further 
works. 

The model developed in this paper can also be directly 
applied to several topics. First, following [H, [g^I, we can 
use our predictions for the 3D density field power spec- 
trum to obtain the power spectrum of the wcak-lensing 
shear (in the Born approximation) , by integrating along 
the line of sight. Second, as in (6ll. |62|| . we can study less 
standard scenarios that involve modifications of gravity 
on cosmological scales, or consider the possible impact 
on the power spectrum of a warm dark matter compo- 
nent [63] . It should also be possible to study the effect of 
neutrinos 64 1 or of baryons [T2j on the power spectrum. 
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Appendix A: Halo model details 

We show in Fig. [^D] the impact on the matter power 
spectrum of some details of the halo model. 

The solid curves show the relative deviation from our 
prediction that would be obtained by using a factor 
{uli - W^), as in instead of {um - W)^ in Eq.(|49l). 
This only gives a k'^ decay at low k instead of a /c^ tail 
(i.e., it satisfies matter conservation but not momentum 
conservation [l^l)- This would slightly overestimate the 
power spectrum on weakly nonlinear scales, where the 
one-halo term starts being non-negligible while remain- 
ing sensitive to its asymptotic low-fc behavior. 

The dashed curves show the relative deviation from 
our prediction that would be obtained by neglecting the 



virial damping factor 



„ikAx\ 



'Aq 



in Eq. 



This would 



overestimate the power spectrum on small scales, in the 
highly nonlinear regime, which are sensitive to motions 
within halo cores. However, on these scales our model 
suffers from other uncertainties, due to the limited accu- 
racy of halo profiles and concentration parameters, and 
to the approximate form of our two-halo contribution. 
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